home *** CD-ROM | disk | FTP | other *** search
/ Collection of Tools & Utilities / Collection of Tools and Utilities.iso / basic / ubas830.zip / DETP.ASM < prev    next >
Assembly Source File  |  1990-01-28  |  6KB  |  382 lines
  1. ;detp.asm
  2.  
  3.     ;determinant of matrix det P
  4.     ;1990 by Yuji KIDA
  5.  
  6. comment &
  7.     usage :
  8.  
  9.     dim Detp%(500):bload "detp",Detp%()
  10.     N=(P-1)\2-2:NN=(N+7)\8*8
  11.     dim A%(N-1,NN-1),W%(N-1,NN-1) ;column size = multiple of 8
  12.     call Detp%(W%(0,0),N,PP,DP)
  13.  
  14.     original BASIC program is :
  15.  
  16.   640   *DETP(P):'           assume every elements >=0 and <p
  17.   650      local I%,J%,K%,DP,A
  18.   660   DP=1:'                                     determinant mod p
  19.   670   for J%=0 to N-2
  20.   680      for I%=J% to N-1:'                      search non_zero
  21.   690         if W(I%,J%) then cancel for:goto 720
  22.   700      next
  23.   710      cancel for:DP=0:goto 840
  24.   720      if I%<>J% then neg DP
  25.   730        :swap block W(I%,J%..N-1),W(J%,J%..N-1)
  26.   740      A=W(J%,J%)
  27.   750      DP=DP*A@P
  28.   760      Inv=modinv(W(J%,J%),P)
  29.   770      for I%=J%+1 to N-1
  30.   780         C=W(I%,J%)
  31.   790         if C then C=C*Inv@P
  32.   800           :for K%=J%+1 to N-1:W(I%,K%)=(W(I%,K%)-C*W(J%,K%))@P:next
  33.   810      next
  34.   820   next
  35.   830   DP=DP*W(N-1,N-1)@P
  36.   840   return(DP)
  37. &
  38.     
  39.     INCLUDE    UB.MAC
  40.  
  41.  
  42.     jmp    start0
  43.  
  44.     even
  45. base_segment    dw    ?
  46. matrixsize    dw    ?
  47. modulus        dw    ?
  48. segstep        dw    ?
  49. var_J        dw    ?
  50. invword        dw    ?
  51. multiplier    dw    ?
  52. counter        dw    ?
  53. detoffset    dw    ?
  54. detsegment    dw    ?
  55. det        dw    ?
  56.  
  57.  
  58.  
  59. modinvword:
  60. ; inp : ax (GCD with 'modulus' must be 1)
  61. ; out : ax = 1/ax mod 'modulus'
  62. ; destroy : nothing
  63.  
  64. ;    cmp    ax,1
  65. ;    je    modinvwordret
  66.     push    bx
  67.     push    dx
  68.     push    si
  69.     push    di
  70.  
  71.     cmp    ax,bx        ;bx=modulus
  72.     jb    modinvw10
  73.     xor    dx,dx
  74.     div    bx
  75.     mov    ax,dx
  76.     cmp    ax,1
  77.     je    modinvwout
  78. modinvw10:
  79.     xchg    ax,bx        ;##
  80.     xor    si,si        ;coef1
  81.     mov    di,1        ;coef2
  82. modinvw20:
  83.     xor    dx,dx
  84.     div    bx
  85.     push    dx        ;remainder
  86.     push    di        ;coef2
  87.     mul    di
  88.     div    cs:[modulus]
  89.     mov    ax,si
  90.     sub    ax,dx
  91.     jae    modinvw30
  92.     add    ax,cs:[modulus]
  93. modinvw30:
  94.     mov    di,ax        ;new coef2=coef1-quotient*coef2
  95.     pop    si        ;new coef1=old coef2
  96.     mov    ax,bx
  97.     pop    bx
  98.     cmp    bx,1
  99.     jne    modinvw20    ;GCD must 1 otherwise endlessloop
  100.     mov    ax,di
  101. modinvwout:
  102.     pop    di
  103.     pop    si
  104.     pop    dx
  105.     pop    bx
  106. modinvwordret:
  107.     ret
  108.  
  109.  
  110. getvalue:
  111.     lodsw
  112.     cmp    ax,1
  113.     jne    getvalueilg
  114.     lodsw
  115.     clc
  116.     ret
  117.  
  118. getvalueilg:
  119.     stc
  120.     ret
  121.  
  122.  
  123.  
  124. detpilg:
  125.     mov    bx,AR1
  126.     mov    word ptr cs:[bx],8001h    ;error
  127.     jmp    detpret
  128.  
  129.  
  130.  
  131. start0:
  132.     mov    ax,cs
  133.     mov    ds,ax
  134.     mov    es,ax
  135.  
  136.     lds_si    v2
  137.     call    getvalue
  138.     jc    detpilg
  139.     mov    cs:[matrixsize],ax
  140.  
  141.     lds_si    v3
  142.     call    getvalue
  143.     jc    detpilg
  144.     mov    cs:[modulus],ax
  145.  
  146.     lds_si    v4
  147.     mov    cs:[detoffset],si
  148.     mov    cs:[detsegment],ds
  149.  
  150.     mov    ax,cs
  151.     mov    ds,ax
  152.     mov    es,ax
  153.  
  154.     mov    [det],1        ;set initial det
  155.  
  156.     mov    bx,V1+2
  157.     mov    ax,[bx]        ;segment of W()
  158.     mov    ds,ax
  159.     xor    si,si
  160.     mov    ax,[si]        ;ax=dimension
  161.     cmp    ax,2
  162.     jne    detpilg
  163.     mov    si,8
  164.     mov    ax,[si]        ;size of 2 dim -1
  165.     inc    ax
  166.     xor    dx,dx
  167.     mov    cx,8
  168.     div    cx
  169.     or    dx,dx
  170.     jnz    detpilg
  171.     mov    cs:[segstep],ax
  172.  
  173.     mov    ax,ds
  174.     add    ax,arrayheadseg
  175.     mov    ds,ax
  176.     mov    es,ax
  177.  
  178.     mov    cs:[base_segment],ax
  179.  
  180.     mov    bx,cs:[modulus]
  181.     mov    cx,cs:[matrixsize]
  182. detp10:
  183.     push    cx
  184.     xor    si,si
  185.     mov    cx,cs:[matrixsize]
  186. detp20:
  187.     mov    ax,[si]
  188.     test    ah,80h
  189.     jz    detp30
  190.     and    ah,7fh
  191.     neg    ax
  192. detp30:
  193.     add    ax,7fffh
  194.     xor    dx,dx
  195.     div    bx    ;modulus
  196.     mov    [si],dx
  197.     add    si,2
  198.     loop    detp20
  199.  
  200.     mov    ax,ds
  201.     add    ax,cs:[segstep]
  202.     mov    ds,ax
  203.     mov    es,ax
  204.  
  205.     pop    cx
  206.     loop    detp10
  207.  
  208.     ;  670   for J%=0 to N-2
  209.  
  210.     mov    ax,cs:[base_segment]
  211.     mov    ds,ax
  212.     xor    si,si        ;ds:si = adr of W(0,J)
  213.  
  214.     mov    ax,cs:[matrixsize]
  215.     dec    ax
  216.     jnz    detp105
  217.     jmp    detp200
  218. detp105:
  219.     mov    cs:[counter],ax
  220.     mov    cs:[var_J],0
  221.  
  222.     ;main loop
  223. detp110:
  224.     mov    ax,ds
  225.     mov    es,ax
  226.  
  227.     ;  680      for I%=J% to N-1:'                      search non_zero
  228.     ;  690         if W(I%,J%) then cancel for:goto 720
  229.     ;  700      next
  230.     ;  710      cancel for:DP=0:goto 840
  231.  
  232.     mov    cx,cs:[matrixsize]
  233.     sub    cx,cs:[var_J]
  234.     xor    ax,ax
  235. detp120:
  236.     cmp    es:[si],ax
  237.     jne    detp130
  238.     mov    dx,es
  239.     add    dx,cs:[segstep]
  240.     mov    es,dx
  241.     loop    detp120
  242.     jmp    detpequ0    ;det=0
  243.  
  244.     ;  720      if I%<>J% then neg DP
  245.  
  246. detp130:
  247.     mov    ax,ds
  248.     mov    dx,es
  249.     cmp    ax,dx
  250.     je    detp150
  251.  
  252.     mov    ax,bx        ;neg det
  253.     sub    ax,cs:[det]    ;
  254.     mov    cs:[det],ax    ;
  255.  
  256.     ;  730        :swap block W(I%,J%..N-1),W(J%,J%..N-1)
  257.  
  258.     mov    cx,cs:[matrixsize]    ;exchange rows
  259.     sub    cx,cs:[var_J]
  260.     mov    di,si
  261. detp140:
  262.     mov    ax,es:[di]
  263.     xchg    [di],ax
  264.     mov    es:[di],ax
  265.     add    di,2
  266.     loop    detp140
  267.  
  268.     mov    ax,ds
  269.     mov    es,ax
  270.  
  271.     ;  740      A=W(J%,J%)
  272.     ;  750      DP=DP*A@P
  273. detp150:
  274.     mov    ax,[si]        ;calc new det
  275.     mul    cs:[det]
  276.     div    bx        ;##
  277.     mov    cs:[det],dx
  278.  
  279.     ;  770      for I%=J%+1 to N-1
  280.  
  281.     mov    cx,cs:[matrixsize]
  282.     sub    cx,cs:[var_J]
  283.     dec    cx
  284.     jcxz    detp190
  285.  
  286.     ;  760      Inv=modinv(W(J%,J%),P)
  287.  
  288.     mov    ax,[si]
  289.     cmp    ax,1
  290.     je    detp155
  291.     call    modinvword
  292. detp155:
  293.     mov    cs:[invword],ax
  294.  
  295.     ;  780         C=W(I%,J%)
  296.     ;  790         if C then C=C*Inv@P
  297.  
  298. detp160:
  299.     mov    ax,es
  300.     add    ax,cs:[segstep]
  301.     mov    es,ax
  302.  
  303.     ;740
  304.  
  305.     mov    ax,es:[si]
  306.     or    ax,ax
  307.     jz    detp185
  308.  
  309.     ;  800           :for K%=J%+1 to N-1:W(I%,K%)=(W(I%,K%)-C*W(J%,K%))@P
  310.     ;                :next
  311.     ;  810      next
  312.  
  313.     push    cx
  314.     mul    cs:[invword]
  315.     div    bx        ;##
  316.     mov    ax,dx
  317.     or    ax,ax
  318.     jz    detp170
  319.     neg    ax
  320.     add    ax,bx        ;##
  321. detp170:
  322.     mov    cs:[multiplier],ax
  323.  
  324.     mov    cx,cs:[matrixsize]
  325.     sub    cx,cs:[var_J]
  326.     mov    di,si
  327. detp180:
  328.     mov    ax,[di]
  329.     mul    cs:[multiplier]
  330.     div    bx        ;##
  331.  
  332.     mov    ax,dx
  333.     xor    dx,dx
  334.     add    ax,es:[di]
  335.     adc    dx,0
  336.     div    bx        ;##
  337.     mov    es:[di],dx
  338.     add    di,2
  339.     loop    detp180
  340.  
  341.     pop    cx
  342. detp185:
  343.     loop    detp160
  344.  
  345.     ;  820   next
  346.  
  347. detp190:    
  348.     inc    cs:[var_J]
  349.     mov    ax,ds
  350.     add    ax,cs:[segstep]
  351.     mov    ds,ax
  352.     add    si,2
  353.  
  354.     dec    cs:[counter]
  355.     jz    detp200
  356.     jmp    detp110
  357.  
  358. detpequ0:
  359.     mov    cs:[det],0
  360.  
  361.     ;  830   DP=DP*W(N-1,N-1)@P
  362.  
  363. detp200:
  364.     mov    ax,[si]        ;final element
  365.     mul    cs:[det]
  366.     div    bx
  367.     mov    ax,dx
  368.     lds    si,dword ptr cs:[detoffset]
  369.     or    ax,ax
  370.     jz    detp210
  371.     mov    word ptr [si],1
  372.     add    si,2
  373. detp210:
  374.     mov    [si],ax
  375.  
  376.     ;  840   return(DP)
  377.  
  378. detpret:
  379.     return
  380.  
  381.  
  382.